1. Field of the Invention
This invention relates to the in-vivo determination and display of estimates of the cardiac ejection fraction, or the end diastolic volume, or both.
2. Description of the Related Art
Information about the output of a patient""s heart is very valuable to a surgical team operating on the patient or to physicians who are trying to diagnose an illness or monitor the patient""s condition. Few hospitals are therefore without some form of conventional equipment to monitor cardiac output.
One common way to determine cardiac output is to mount some flow-measuring devices on a catheter, and then to thread the catheter into the patient and to maneuver it so that the devices are in or near the patient""s heart. Some such devices inject either a bolus or heat at an upstream position, such as in the right atrium, and determine flow based on the characteristics of the injected material or energy at a downstream position, such as in the pulmonary artery.
For example, U.S. Pat. No. 4,236,527 (Newbower et al., Dec. 2, 1980) and U.S. Pat. No. 4,507,974 (Yelderman, Apr. 2, 1985), describe systems for measuring cardiac output in which heat is used as an indicator. In such heat-based systems, a balloon catheter is typically positioned proximal to the branch of the pulmonary artery via the right atrium and the right ventricle. The catheter includes a resistive heating element, which is positioned in the atrium and/or ventricle, and a thermistor, which is positioned in the artery. Cardiac output is then calculated as a function of the sensed downstream temperature profile.
U.S. Pat. No. 5,146,414 (McKown, et al., Sep. 8, 1992) describes a system in which the transfer function of the channel (the region from where an indicator such as heat is applied to the blood upstream to the downstream position where the indicator concentration, such as temperature, is sensed) is modeled, the approximate spectrum of the noise is determined, and the output of the system is used in a feedback loop to adaptively update the parameters of the model and thus to improve the estimate of cardiac output (CO). U.S. Pat. No. 5,687,733 (McKown, et al., Nov. 18, 1997) describes an improvement over the earlier McKown ""414 system that estimates both the CO trend and an instantaneous CO value. Moreover, in the McKown systems, only the zero-frequency (dc or steady state) gain of the channel is required to get an estimate of the cardiac output (CO).
Although these known systems provide estimates of cardiac output with varying degrees of accuracy, they fail to provide any estimate of the heart""s ejection fraction (EF), typically, the right ejection fraction (REF), which is defined as the ratio between the stroke volume (SV) of the heart and its end diastolic volume (EDV). The ejection fraction is thus a measure of how efficiently the heart pumps out the blood that it can contain.
Because of its diagnostic importance, there are several known methods for measuring EF. Such systems, however, frequently rely on the use of an injected bolus and on evaluation of the washout (thermodilution) curve in the blood vessel. U.S. Pat. No. 4,858,618 (Konno, et al., issued Aug. 22, 1989), for example, describes a thermodilution system for determining right ventricular ejection fraction. In this known system, a cold bolus indicator is injected into the right ventricle. Pre- and post-bolus temperatures are sensed in the pulmonary artery. The temperature differentials are used to determine the ejection fraction.
One problem with using a bolus to determine EF is that it is difficult to establish just where on the sensed bolus curve the measurements are to begin, since the front side of the curve depends heavily on mixing, on the heart rate, and even on how fast the administering nurse is pushing the syringe plunger while injecting the bolus. Another problem faced by all such known systems is that they require synchronization with the heart cycle in order to reduce the effects of the heartbeat when producing an EF estimate. Some systems synchronize based on plateaus in the washout curve, but this presupposes a fast and very accurate thermistor. Other systems rely for synchronization on an EKG trigger. EKG synchronization, however, is difficult, since it is then necessary to slave in and precisely coordinate the timing of other instruments, each gathering its own data.
Further problems of existing systems for determining EF stem from their need to identify discrete plateaus in the dilution profiles created by the heart beats. This is necessary because these systems use the plateaus as markers in order to fit exponential or ratio-based curves to the data, which are in turn used to evaluate the dilution decay. This approach is accurate in practice, however, only for a relatively slow heart rate and a thermistor whose response is significantly faster than the decay parameter xcfx84.
In effect, these conventional systems assume a square-wave dilution curve. This is, however, usually an unrealistic assumption. First, most of the patients needing EF measurements in a hospital are not in the best of health; rather, they tend to have relatively high and erratic heart rates. Furthermore, in systems that use a bolus of relatively cold fluid, the sensed heart rate is likely to be incorrect since the cold bolus itself tends to affect not only the heart rate, but also its regularity. Second, real thermistors distort the plateaus, so that the exponential fits themselves become distorted. Third, as the EF rises, the drops in the plateaus also rise. This causes the systems to use fewer plateaus, and thus reduce their accuracy, because of the limited signal-to-noise ratios of these systems.
For example, one known system uses a fast response injectate cardiac output pulmonary artery catheter together with an electrocardiogram R-wave detector to measure EF and EDV. The exponential method of measuring REF then synchronizes R-wave events with plateaus occurring during the downslope of a thermodilution curve and fits the decay of the curve with an exponential function. Thus, if T(i) is the PA temperature after the i-th R-wave and T(ixe2x88x92n) is the temperature n R-waves earlier in time, then:
T(i)=T(ixe2x88x92n)*exp(xe2x88x92t/xcfx84),xe2x80x83xe2x80x83(Equation 1)
where t is time and xcfx84 is the decay parameter.
The physiological washout decay can then be represented by (1xe2x88x92EF)n, where n is the number of R-waves in the observation interval (for example, from 80% down to 30% of the peak). One can then represent time in terms of the heart rate (HR):
t=n*60/HRxe2x80x83xe2x80x83(Equation 2)
where HR is the local average from the (ixe2x88x92n)""th to the i""th R-wave in beats per minute. Given these relationships, the following can then be shown:
EF=1xe2x88x92exp(xe2x88x9260/(xcfx84*HR)).xe2x80x83xe2x80x83(Equation 3)
One of the problems with this system is that the thermistor must have a sufficiently fast response time to allow measurement of the true physiological decay time. At low heart rates, this puts plateaus in the temperature data during systole, which must be dealt with in determining the decay parameter T. This is, indeed, the primary reason for the R-wave synchronization, since, other than that, the local average HR is all that is required.
Another problem with this known system is that it is bolus-based and intermittent by nature. In addition, only part of the temperature data is used (from the R-wave around 80% washout to the R-wave around 30% washout: typically one-five R-waves). This introduces variability or lack of precision into the measurement of the injectate cardiac output (ICO) due to irregular R-wave intervals or large noise sources such as respiratory ventilation.
The earlier McKown systems improve on such a bolus-based approach by instead generating an input injectate signal, preferably in the form of a pseudo-random binary heat signal, and then estimating the parameters of a transfer function model of the input-output channel. The preferred model used is the lagged normal transfer function (described below). Both the measurement and the modeling of the transfer function model are carried out in the frequency domain at the harmonics of the preferred input injectate signal. In order to understand the weaknesses of these systems, it is helpful to have at least a basic understanding of the lagged normal model of the transfer function.
In the context of estimating cardiac output, the xe2x80x9clagged normal modelxe2x80x9d described by Bassingthwaighte, et al. in xe2x80x9cApplication of Lagged Normal Density Curve as a Model for Arterial Dilution Curves,xe2x80x9d Circulation Research, vol. 18, 1966, has proven to be particularly accurate and useful, and it is therefore the model for cardiac output used in, for example, McKown ""733. The lagged normal model is defined as a linear, time-invariant system (LTIS) whose impulse response is the convolution of a unity-area Gaussian (normal distribution) function and a unity-area decaying exponential. The Gaussian has two parameters: the mean xcexc and the standard deviation "sgr". The exponential has one parameter: the time-decay parameter xcfx84. The unity-gain, lagged-normal transfer function Hxe2x80x94LN at each frequency xcfx89 sampled (xcfx89 is the independent variable in this model) thus depends on xcexc, "sgr", and xcfx84 as follows:
Hxe2x80x94LN(xcfx89|xcexc,"sgr",xcfx84)=exp[xe2x88x92jxc2x7xcfx89xc2x7xcexcxe2x88x92(xcfx89xc2x7"sgr")2/2]/(1+jxc2x7xcfx89xc2x7xcfx84)xe2x80x83xe2x80x83(Equation 4)
where exp is the exponential function and the physical meaning of the parameters is:
xcexc: a pure time delay that represents translational flow
"sgr": a measure of random dispersion
xcfx84: the decay parameter, that is, a time constant associated with mixing in a distribution volume, which, in this example, is the blood vessel.
The units of xcexc, "sgr", and xcfx84 are time (seconds) and the units of xcfx89 are radians per second.
This model is used in, for example, not only in the McKown ""733 system, but also in a more recent system described in the U.S. patent application Ser. No. 09/094,390, FILED ON Jun. 9, 1998 by the same inventors as those of the present invention, which build on the McKown ""733 techniques.
Although other indicators may be used, in the preferred embodiment of these systems, heat is used as the indicator, and the indicator driver signal is a pseudo-random binary sequence (PRBS). The driver/sensor pair therefore preferably consists of a heater and a thermistor. Hxe2x80x94LN is then estimated as an optimized fitting of a vector of complex values Hxy(xcfx89n), each representing a measurement of the transfer function between a heater power signal x and a thermistor temperature signal y. Each vector contains the parameters fitted to the measured temperature data at each of ten frequencies xcfx89n (the first ten PRBS harmonics).
More specifically, the system in McKown ""733 computes the state vector X=[dc, xcexc, "sgr", xcfx84] that minimizes the cost function Cost_Hxy, defined as:
Cost_Hxy=SUM[Hxyxe2x80x94SAE(xcfx89n|X)xc2x7W(xcfx89N)]xe2x80x83xe2x80x83(Equation 5)
where the sum is taken over N=1 to 10 (or however many harmonics are used), W(xcfx89n) are weights and:
Hxyxe2x80x94SAE(xcfx89n|X)=[Hxy_avg(xcfx89n)xe2x88x92dcxc2x7Hxyxe2x80x94LN(xcfx89N|X)]2xe2x80x83xe2x80x83(Equation 6)
which is the squared absolute error (SAE) of the averaged measured transfer function Hxy_avg(xcfx89n) relative to the lagged normal transfer function model Hxyxe2x80x94LN(xcfx89n|X) at the PRBS harmonic frequencies xcfx89n given the state vector X=[dc, xcexc, "sgr", xcfx84].
Once xcexc, "sgr", and xcfx84 are known, then each of the ten complex measured numbers Hxy(xcfx89N) would individually provide an estimate of cardiac output (CO) according to:
CO(n)=Kxc2x7Hxe2x80x94LN(xcfx89n)/Hxy(xcfx89n) for n=1 to 10xe2x80x83xe2x80x83(Equation 7)
where K is a known or experimentally determinable conversion constant.
In order to apply this relationship, the McKown ""733 system first determines not only what the values of xcexc, "sgr", and xcfx84 should be, but also how the ten cardiac output estimates CO(n) should be combined. One should note that the cardiac output does not depend on the shape of H(xcfx89), or Hxy(xcfx89), but only on the zero-frequency gain, dc, of Hxy. Since the experimental transfer function Hxy is measured at ten frequencies xcfx89n that are not zero, however, the McKown ""733 system in essence extrapolates the measured Hxy(xcfx89) to zero frequency. A known optimization routine is then used to provide a best fit of the ten modeled transfer function values H_xy to the observed values. The relationship shown above for CO can then be reduced to CO=K/dc, where dc is the zero-frequency (xcfx89=0) gain value in units of degrees Celsius per watt, and K is an experimentally determined constant with the unit (liters per minute)/(degrees Celsius per watt).
Note that the McKown ""733 system provides a continuous CO value (equivalently, the dc value), as well as the decay parameter xcfx84. Note that xe2x80x9ccontinuousxe2x80x9d does not here mean that displayed values are xe2x80x9ccontinuously changing,xe2x80x9d but rather that they can be updated every processing cycle (preferably a PRBS cycle), after an initialization period.
There are, however, problems with the prior-art technologies, which are based on frequency domain (typically, cross-correlation) transfer function measurement and modeling. A primary limitation of these prior-art techniques is that stated by W. D. Davies in xe2x80x9cSystem Identification for Self-Adaptive Control,xe2x80x9d Wiley-Interscience, 1970, namely, that xe2x80x9csince the technique described here may also be considered as one of identifying the frequency response of an unknown system, it will also unfortunately combine in the final estimate the frequency components of the noise that lie within the system bandwidth, and to date there exists no theory that allows the separation of the signal from the noise.xe2x80x9d
In the McKown ""414 and ""733 systems, for example, only the transfer function""s dc-gain is used, whereas, in other systems, such as in the U.S. patent application Ser. No. 09/094,390, the indicator decay time-constant xcfx84 is used, in addition to a heart rate estimate. A problem of these prior art approaches is, however, the degree to which estimation errors are coupled between the parameters dc, xcfx84, "sgr", and xcexc. This coupling is primarily due to the low-frequency thermal (indicator) noise, for example, the noise created by the patient""s respiration, whether natural or mechanically ventilated. The parameter estimation then adversely affects the estimate of the cardiac output, and also often degrades the accuracy of the estimated REF and EDV so badly that the measurements become clinically unacceptable.
Another problem of these earlier techniques is their use of the four-parameter (dc, xcexc, "sgr", xcfx84) lagged normal frequency domain model to analyze the transfer function data. Typically, if enough noise is present, then the optimization routine (for example, squared error cost function minimization) may converge to local or false minima for the vector (xcexc, "sgr", xcfx84) of shape parameters. In other words, there may be several xe2x80x9cbestxe2x80x9d combinations of xcexc, "sgr", and xcfx84, most or even all of which are bad in the sense of lying too far from the true values. Although this affects the quality of continuous cardiac output (CCO) measurements only slightly (due to some dcxe2x88x92xcfx84 coupling) it is a major hindrance to the accurate determination by existing systems of continuous EF/EDV since an accurate estimate of xcfx84 is required.
One other shortcoming of the frequency-domain techniques described above that use the lagged normal model is that they calculate estimates based on only a limited number of harmonics. Consequently, faster time constants, which lie outside the bandwidth of the primary (first ten) PRBS harmonics, are poorly determined.
What is needed is therefore a system that can produce continuous estimates of the EF or EDV, or both, but whose estimates are substantially unaffected by low frequency-induced errors in the dc and xcfx84 estimates; this in turn would provide more accurate CO and EF/EDV measurements, respectively. This invention provides such a system, and a corresponding method, for determining CO and EF/EDV.
The invention provides a method for estimating a cardiac performance value, such as cardiac output CO and/or cardiac ejection fraction EF of a patient according to which an indicator (preferably heat) is injected at an upstream position in a patient""s heart according to a predetermined injected indicator signal x(t). An indicator concentration sensor such as a thermistor then senses a local indicator concentration signal y(t) at a downstream position. The region from and including the upstream position to and including the downstream position forms a channel for the blood.
The indicator concentration signal is then divided into at least one sub-signal that is synchronous with the injected indicator signal. A processing system then calculates a first time-domain, channel relaxation model that has each sub-signal as an input. It then also calculates a decay parameter xcfx84 as a pre-determined function of the first, time-domain channel relaxation model. The processor then estimates the cardiac performance value as a predetermined function of the decay parameter xcfx84.
In the preferred embodiments of the invention, the injected indicator signal x(t) is generated as a series of alternating transitions between a high state and a low state. Examples of suitable injected indicator signals includes a periodic, pseudo-random binary sequence (the preferred embodiment), trains of random or periodic random square waves, and even non-binary signals such as trigonometric functions and spread spectrum signals.
According to one aspect of the invention referred to as y_tau integration, for each transition of the injected indicator signal, a corresponding segment of the indicator concentration signal is isolated, with each segment comprising one of the sub-signals. For each segment of the indicator concentration signal, a segment relaxation parameter is then calculated. The processor then calculates the decay parameter xcfx84 as a predetermined function of the segment relaxation parameters.
In implementations of the invention in which the injected indicator signal is periodic, with a plurality of transitions during each period, each sub-signal of the indicator concentration signal corresponds to one period of the injected indicator signal. The y_tau integration embodiment of the invention, further preferably includes sign-rectification of all the segments before the decay parameter xcfx84 is calculated.
In y_tau integration, the step of calculating the decay parameter xcfx84 preferably includes the sub-steps of generating the first time-domain, channel relaxation model as a time-domain exponential function of the decay parameter; calculating a cost function that is a predetermined function of the sum of differences between the exponential function of the decay parameter and the respective segments of the indicator concentration signal; and calculating the decay parameter xcfx84 by determining a minimum of the cost function.
In embodiments of the invention in which CO is to be estimated, the system according to the invention includes a heart rate monitor. The cost function is then preferably a predetermined function of both the decay parameter xcfx84 and a steady-state channel gain parameter (dc). The processing system according to the invention then determines optimum values of the decay parameter xcfx84 and the steady-state channel gain parameter dc that minimize the cost function. It then calculates the cardiac output (CO) value as a predetermined function of the optimum value of the steady-state channel gain parameter, and calculates the cardiac ejection fraction (EF) as a predetermined function of the optimum value of the steady-state channel gain parameter and the measured heart rate (HR).
According to a second aspect of the invention referred to as y_avg integration, the processor divides the indicator concentration signal y(t) into a plurality of the sub-signals such that each sub-signal corresponds to one period of the injected indicator signal. An average of the sub-signals is then calculated to form an averaged indicator concentration signal. The channel relaxation model is thereby generated as a time-domain, lagged-normal function of both the decay parameter xcfx84 and the steady-state channel gain parameter (dc). A cost function is then evaluated that is a predetermined function of the difference between the averaged indicator concentration signal and the time-domain, lagged-normal function convolved with the injected indicator signal. The system then determines optimum values of the decay parameter xcfx84 and the steady-state channel gain parameter dc that minimize the cost function. Values for CO and EF are then calculated as predetermined functions of the optimum values of the steady-state channel gain parameter (for CO) and the steady-state channel gain parameter and the measured heart rate (for EF).
In another embodiment of the invention that includes combined parameter estimation, estimates of the decay parameter xcfx84 and of the steady-state gain dc are determined in three different ways: using y_tau integration, y_avg integration, and by determining the optimum values to minimize a cost function based on a frequency-domain lagged normal model of the channel. The three different estimates are then normalized and combined using a weighted average.